A strong nonrandom pattern in Matlab default random number generator

The default random number generator in Matlab versions between 5 and at least 7.3 (R2006b) has a strong dependence between the numbers zi+1, zi+16, zi+28 in the generated sequence. In particular, there is no index i such that the inequalities zi+1 < 1/4, 1/4 ≤ zi+16 < 1/2, and 1/2 ≤ zi+28 are satisfied simultaneously. This fact is proved as a consequence of the recurrence relation defining the generator. A random sequence satisfies the inequalities with probability 1/32. Another example demonstrating the dependence is a simple function f with values −1 and 1, such that the correlation between f(zi+1, zi+16) and sign(zi+28 − 1/2) is at least 0.416, while it should be zero. A simple distribution on three variables that closely approximates the joint distribution of zi+1, zi+16, zi+28 is described. The region of zero density in the approximating distribution has volume 4/21 in the three dimensional unit cube. For every integer 1 ≤ k ≤ 10, there is a parallelepiped with edges 1/2, 1/2 and 1/2, where the density of the distribution is 2. Numerical simulation confirms that the distribution of the original generator matches the approximation within small random error corresponding to the sample size.